function logpibpZ = logpibp(Z, alphav)
Z = Z ~= 0;
Z = Z(:,any(Z));

[N,K] = size(Z);
logpibpZ = K*log(alphav) - K*gammaln(N+1) - alphav*hrmsum(N);
feahist = containers.Map('KeyType', 'char', 'ValueType', 'double');
for k = 1:K
    histr = sprintf('%d', Z(:,k));
    if feahist.isKey(histr)
        feahist(histr) = feahist(histr) + 1;
    else
        feahist(histr) = 1;
    end
    mk = sum(Z(:,k));
    logpibpZ = logpibpZ + gammaln(N-mk+1) + gammaln(mk);
end
logpibpZ = logpibpZ - sum(cellfun(@(x)gammaln(x+1), feahist.values));

end